author=Marcio name=Listas ===== Fisica Computacional ===== Professor: [[http://complex.if.uff.br/marcio|Marcio Argollo]] Contato: Sala A1-02, tel. 2629-5773 **Referências básicas:** * Manuais de Fortran 77: * http://www.stanford.edu/class/me200c/tutorial_77/ * http://www.math.hawaii.edu/~hile/fortran/fortmain.htm * http://www.fortran.com/F77_std/rjcnf0001.html * http://www.idris.fr/data/cours/lang/fortran/f90/F77.html * Manuais de Fortran 90 * http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_fortran90.pdf * http://www.cisl.ucar.edu/tcg/consweb/Fortran90/F90Tutorial/tutorial.html * Fortran 90 para programadores de Fortran 77 * http://www.nsc.liu.se/~boein/f77to90/ * Manuais de C/C++ * http://imada.sdu.dk/~svalle/courses/dm14-2005/mirror/c/ * http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_C.pdf * http://www.cs.cf.ac.uk/Dave/C/CE.html * http://www.ic.unicamp.br/~cmrubira/aacesta/cpp/cpp15.html * Livros-texto * [[http://sip.clarku.edu/|An Introduction to Computer Simulation Methods: Applications to Physical Systems]], H. Gould e J. Tobochnik. * Física em computadores, Paulo Murilo Castro de Oliveira e Suzana Maria Moss de Oliveira, São Paulo: Livraria da Física, 2000. (Cópia na biblioteca) * Números no computador: * Representação binária: http://en.wikipedia.org/wiki/Binary_numeral_system * Números positivos e negativos em binário: http://www.cs.cornell.edu/~tomf/notes/cps104/twoscomp.html * Números reais e representação de ponto-flutuante: http://www.psc.edu/general/software/packages/ieee/ieee.php e http://www.cs.cornell.edu/~tomf/notes/cps104/floating.html O aluno é fortemente estimulado a procurar fontes de informação online (cursos, vídeo-aulas, etc.) para complementar os livros-texto. Para uma revisão de comandos básicos em Unix/Linux, utilização do gnuplot e programação em C veja o tutorial apresentado na Semana da Física 2009 [[http://complex.if.uff.br/marcio/semana2009|"Ferramentas de Computação em Física"]]. ================================================================================= ==== Listas de Exercícios ==== As listas deverão ser entregues por email, compactadas no formato tar.gz e com o nome listaN_seunome.tar.gz Para compactar seus arquivos (por exemplo arq1.c e arq2.eps no formato tar.gz), gerando o arquivo arquivofinal.tgz, digite marcio@pelego:~/cursos/semana2009$ tar -zcvf arquivofinal.tar.gz arq1.c arq2.eps marcio@pelego:~/cursos/semana2009$ Para compactar todo um diretorio (por exemplo meudir), digite marcio@pelego:~/cursos/semana2009$ tar -zcvf arquivofinal.tar.gz meudir marcio@pelego:~/cursos/semana2009$ === Lista 1) Exercicios de revisão de programação - Data de entrega: 18/08/2011 (Extendido até 24/08/2011) === - Compute os 50 primeiros números de Fibonacci F_n=F_{n-1}+F_{n-2} e faça um gráfico F_n/F_{n-1} versus n, mostrando que a razão converge para o número áureo \phi \approx 1.6180339887. Veja mais sobre Fibonacci e a razão áurea [[http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/|aqui]]. - Somas não são comutativas no computador! Calcule as somas S(N)=\sum_{i=1}^{N-1} \frac{1}{i^2} e SI(N)=\sum_{i=1}^{N-1} \frac{1}{(N-i)^2} com i sendo uma variável ponto flutuante de precisão simples e coloque **em um mesmo gráfico** as diferencas \left|S(N)-S_{\infty}\right|/S_{\infty} vs N e \left|SI(N)-S_{\infty}\right|/S_{\infty} vs N para N=100,200,..,10^6 onde S_{\infty}=\frac{\pi ^2}{6}. - Imagine o empilhamento de dominós de massa unitária e comprimento 2. O primeiro dominó se encontra com o centro de massa na origem, e a cada passo n deslocamos os dominós ja empilhados de \delta_n para a direita e adicionamos outro dominó na parte inferior da pilha, também com seu centro de massa na origem. Escreva um programa que diz, para um dado \delta_n=\delta~~ \forall n, quantos dominós podem ser empilhados até que o centro de massa ultrapasse a largura do primeiro dominó e a pilha tombe. (Veja uma discussão do problema [[http://web.mat.bham.ac.uk/C.J.Sangwin/Teaching/pus/dominoes.pdf|aqui]]). **A partir deste exercício todos os cálculos envolvendo reais devem ser feito em ponto flutuante com precisão dupla** === Lista 2) Derivada numérica e integração de equações diferenciais: O método de Euler - Data de entrega: 06/09/2011 == Sendo a derivada de uma função f(x) definida por \frac{df}{dx}=\lim_{h \to 0}\frac{f(x+h)-f(x)}{h} e f(x+h)=f(x)+hf^{\prime}(x)+\frac{h^2}{2}f^{\prime \prime}(x)+\ldots, o método de Euler (com difereças posteriores) corresponde ao truncamento desta expansão em ordem h. Obtemos assim que a derivada no ponto x é aproximada por \frac{df}{dx}=\frac{f(x+h)-f(x)}{h}. - Mostre que o método de Euler com diferenças centradas é exato em {\cal O}(h^2), onde h é o passo de discretização da derivada. Nesse método escrevemos \frac{df}{dx}=\frac{f(x+h)-f(x-h)}{2h} - Calcule a derivada numérica de f(x)=\frac{sin(x^2)e^{x/3}}{\sqrt{(x^2+4)}} no ponto x=3 usando o método de Euler simples com passos de discretização h=2^{-25},2^{-22},2^{-20},2^{-15},2^{-12},2^{-10},2^{-5},2^{-1} ~e~ 2^{0} e faça um gráfico mostrando a diferença relativa entre a derivada numérica e o valor exato da derivada para variáveis tipo ponto flutuante com precisão simples e dupla. Entre quais valores de h a aproximação é aceitável? (Solução [[http://davinci.if.ufrgs.br/wiki/index.php/Derivada_Num%C3%A9rica|aqui]]) - Mostre que o método de Euler é instável para qualquer passo h na solução do problema de crescimento exponencial \frac{dy}{dx}=\lambda x - Resolva numericamente a equação diferencial \frac{dy}{dt}-y=-\frac{1}{2}e^{t/2}\sin(5t)+5e^{t/2}\cos(5t) com y(0)=1 com o método de Euler simples para diferentes valores de passo h=2^{-20}, 2^{-15}, 2^{-10} e 2^{-5} e compare cada solução nos instantes t=1,2,3,4 e 5 com a solução exata y(t)=e^t+e^{t/2}sin(5t). Voce pode encontrar essa solução exata digitando a equação diferencial no site do [[http://www.wolframalpha.com|WolframAlpha]]. Experimente! - Faça um gráfico comparando a solução numérica com passo h=2^{-5} e a solução exata entre t=0 e t=5. - Resolva a equação acima com o método de Euler com diferenças centradas e os mesmos valores de h. Para obter o primeiro passo use o método de Euler simples com passo de discretização dez vezes menor que o usado no resto do intervalo. * Referências * http://mathworld.wolfram.com/NumericalDifferentiation.html * http://davinci.if.ufrgs.br/wiki/index.php/M%C3%A9todo_de_Euler === Lista 3) Equações diferenciais ordinárias de sistemas conservativos: métodos simpléticos - Data de entrega: 08/09/2011== Ver [[http://jcp.aip.org/resource/1/jcpsa6/v97/i3/p1990_s1?isAuthorized=no| Reversible multiple time scale molecular dynamics]] (pdf [[http://www.columbia.edu/cu/chemistry/groups/berne/papers/jcp_97_1990_1992.pdf|aqui]]) Estes métodos procuram criar esquemas de discretização da evolução Hamiltoniana de um sistema conservativo, de modo que o volume de pontos representativos no espaço de fase seja preservado durante a evolução temporal. Estes esquemas envolvem fatorização do operador Liouvilleano. Veja nesse {{:molecular_dynamics.pdf|excelente texto}} (seção 4 página 12) uma discussão sobre o assunto. - Use o método de Euler-Cromer e o método de Euler com passo h=10^{-1} e 10^{-4} para resolver numericamente o problema do pêndulo simples na aproximação de pequenas oscilações. Ulitize a condição inicial \theta_0=5^o e v_0=0. Faça um gráfico contendo as soluções analítica e numérica da a evolução temporal da posição angular por aproximadamente 10 períodos do movimento. Na aproximação numérica use os métodos de Euler e Euler-Cromer. Faça outro gráfico ilustrando a evolução temporal da energia no mesmo alcance temporal para as 3 soluções. === Lista 4) Equações diferenciais ordinárias: Métodos de Runge-Kutta - Data de entrega: 13/09/2011== A aproximação para a derivada, como sugerida no método de Euler, introduz erros sistemáticos que acarretam na destruição de certos invariantes da dinâmica original, como a energia total no movimento do pêndulo simples. Para reduzir os efeitos da aproximação numérica à derivada, pode-se recorrer a métodos com aproximações de ordem superior (o método de Euler é de primeira ordem no passo de discretização h). Para informações básicas sobre o método: http://mathworld.wolfram.com/Runge-KuttaMethod.html \\ **Atenção** \\ Para implementar o método com sistemas de equações diferenciais: http://www.nsc.liu.se/~boein/f77to90/rk.html - Resolva numericamente o problema do pêndulo físico utilizando o método de Runge-Kutta de quarta ordem com passo h=2^{-15} para reproduzir os resultados do artigo [[http://www.sbfisica.org.br/rbef/indice.php?vol=17&num=1|"Comportamento crítico no pêndulo simples"]] do professor Paulo Murilo. Faça um relatório seguindo os passos deste trabalho. (veja {{:pendulo_pmco.tgz|aqui}} um programa em C e outro em Fortran77 para calcular a evolução temporal do pêndulo) === Lista 5: Caos em sistemas dissipativos: o pêndulo forçado-amortecido -Data de entrega:20/09/2011 === - Leiam {{:universality_in_chaos.pdf|essa excelente referência}} sobre caos em sistemas dissipativos do Predrag Cvitanović. Nesta sequência de problemas utilize o método de Runge-Kutta de quarta ordem com passo h=10^{-4} Considere o pêndulo simples em um meio viscoso, que oferece resistência ao seu movimento com uma força proporcional à sua velocidade F_v=-\gamma v. A equação de movimento pode ser escrita como \frac{d^2\theta}{dt^2} + \frac{1}{q}\frac{d\theta}{dt} + \sin{(\theta)}=0.\\ - Faça o gráfico da posição como função do tempo para q=0.2,2.0 e 20.0 com posição angular inicial \theta_0=\pi/2 e \omega_0=0. Note que para valores grandes de q devemos recuperar o movimento sem dissipação (conservativo). - Considere agora uma força periódica impulsionando o pêndulo com força F_e = A\cos{(\omega_e t)}. A equação de movimento pode ser escrita como \frac{d^2\theta}{dt^2} + \frac{1}{q}\frac{d\theta}{dt} + \sin{(\theta)}=A\cos{(\omega_E t)}. - Construa o gráfico \omega(t) \times t com q=2, ~ \omega_E = 2/3 e A=1.08, 1.10, 1.13 e 1.23, respectivamente. Observe a mudança de comportamento ao variarmos a amplitude A da força externa. Simule o pêndulo por 100 ciclos da força externa e despreze os primeiros 10 ciclos (quanto vale o período \tau da força externa neste exemplo?). - Faça um gráfico "estroboscópico", chamado seção de Poincaré, imprimindo os pontos no espaço de fase correspondentes a intervalos de tempo múltiplos do período \tau da força externa t=0,\tau,2\tau,\ldots,n\tau,\ldots. - Obtenha uma série de pontos \omega_0,\omega_1,\ldots,\omega_n,\ldots onde \omega_n corresponde ao valor da velocidade angular quando \theta=0 pela n-ésima vez durante a evolução temporal do pêndulo. Despreze os valores obtidos para tempos menores que 10 ciclos da força externa. Construa com esses pontos outra seção de Poincaré em um gráfico \omega_n \times \omega_{n+1} para os mesmos parâmetros do problema anterior. === Lista 6: Caos no mapa logístico - 27/09/2011 === - Para ler mais sobre o problema: http://mathworld.wolfram.com/LogisticMap.html. - Veja {{:swinney-chaos_experiments.pdf|aqui}} uma referência clássica onde se obtém o mapa logístico em um sistema experimental. E {{:may.pdf|aqui}} uma sobre seu uso no contexto ecologico. Considere o mapa logístico x^{\prime}=f(x)=\lambda x(1-x) - Faça os gráficos (escala lin-log) da evolução temporal de x a partir de um valor inicial arbitrário para \lambda=0.01 e \lambda=0.9. - Fazer "gráficos de escada" para os mesmos valores de parâmetro e condição incial x=0.6. Como critério de parada use a condição x-x^{\prime}<10^{-3}. - Fazer "gráficos de escada" da primeira iterada para o caso \lambda=3.1 e condição inicial x=0.3 e para a segunda iterada f(f(x)) com mesmos \lambda e condição inicial. - Construir o diagrama de bifurcação do mapa logístico na região 0.0<\lambda\le 4.0. - Obter o diagrama de bifurcação para os mapas **a)** f(x)=r\sin{\left(\pi x\right)} **b)** f(x)=\left(27/4\right) r x^2(1-x) ambos na região 0.75 === Lista 7: Caos no mapa logístico II - 06/10/2011 === - Leia sobre a teoria de sistemas dinâmicos {{:chaos_in_dynamical_systems-ott81.pdf|aqui}}. - Artigo de revisão sobre a {{:feigenbaum.pdf|teoria de escala para bifurcações de período}} e caos em mapas (escrito pelo criador da teoria). - Excelente {{:kadanoff83-chaos_in_turbulence.pdf|artigo de revisão}} sobre turbulência e caos. - Obter analiticamente os pontos fixos para órbitas de período 2 e numericamente os de período 4 e analisar as regiões de estabilidade destes pontos fixos. Para encontrar as raízes voces podem usar o [[http://www.wolframalpha.com|WolframAlpha]]. Escreva, por exemplo, //roots(l*x*(1-x)-x)//. Para calcular derivadas escreva, por exemplo //diff(sin(x),x)//. No caso de órbitas de período 4 voces devem escolher o algoritmo numérico para zero de funções de sua preferência. - O mapa logístico com \lambda=4.0 é ergódico, o que significa que podemos substituir médias temporais por médias sobre configurações de um ensemble e, nesse caso, \Gamma(\lambda)=\lim_{n\to \infty}\frac{1}{n}\sum_{i=1}^n \log{\left|f^{\prime}(x_i)\right|}=\int \log{\left|f^{\prime}(x)\right|}\rho(x) dx onde \rho(x) é a fração de vezes que a órbita passa entre x e x+dx. Calcule, subdividindo o intervalo [0,1] em 1000 subintervalos, o histograma \rho(x) para as primeiras 100000 iterações (após um transiente de 1000 iterações) do mapa logístico com \lambda=4.0 e encontre \Gamma(\lambda). Esse caso particular admite a solução analítica \rho(x)=\frac{1}{\pi \sqrt{x(1-x)}}. Compare os resultados mostrando-os em um mesmo gráfico. - Obter os expoentes de Lyapunov \Gamma(\lambda) para os mapas **a)** f(x)=\lambda x(1-x) na região 0.0<\lambda\le 4.0. **b)** f(x)=r\sin{\left(\pi x\right)} **c)** f(x)=\left(27/4\right) r x^2(1-x) ambos b) e c) na região 0.75 === Lista 8: Análise de séries temporais - 11/10/2011 === (1) Voce encontra [[http://cursos.if.uff.br/fiscomp-0210/lib/exe/fetch.php?media=pi_10kdigits.data.gz|neste arquivo]] a sequência dos primeiros 10 mil dígitos de \pi e {{:celegans.tgz|aqui}} a sequência das primeiras 10 mil bases de nucleotídeos compondo o cromossomo I do genoma do verme Caenorhabditis elegans (C.elegans), o primeiro eucarioto a ter seu genoma completamente sequenciado (ver http://users.rcn.com/jkimball.ma.ultranet/BiologyPages/C/Caen.elegans.html e http://en.wikipedia.org/wiki/Caenorhabditis_elegans para mais detalhes). Neste arquivo 0,1,2,3 correspondem a A,T,C e G, respectivamente. (a) Calcular a frequência P(i) de ocorrência de cada símbolo i da sequência de \pi e do genoma. (b) Calcular a frequência P(ij) de ocorrência de cada par de símbolos (ij) da sequência de \pi e do genoma. (2) A função de correlação c(\tau)=\frac{\left(2/T\sum_{t=0}^{T/2}x(t)x(t+\tau)\right ) - \left(2/T\sum_{t=0}^{T/2}x(t)\right )\left(2/T\sum_{t=0}^{T/2}x(t+\tau)\right )}{ \left(2/T\sum_{t=0}^{T/2}x(t)^2\right ) - \left(2/T\sum_{t=0}^{T/2}x(t)\right )^2} é uma das ferramentas tradicionais de análise de sinais. Ela mede o grau de correlação entre o valor do sinal medido em um instante e outro sinal medido em um instante de tempo subsequente. Note que, em sinais não estacionários, a função de correlação pode depender do instante inicial da medida c(t,\tau). Não trataremos desses casos aqui. (a) Calcule a função de correlação para cada uma das séries de 100000 pontos obtidas com a evolução do mapa logístico com \lambda=2.1 (ponto fixo), \lambda=3.0 (ponto de bifurcação), \lambda=3.2 (movimento periódico ou ciclo limite), \lambda=3.569945672 (ponto de acumulação, infinitas órbitas periódicas instáveis emergem neste ponto), \lambda=3.7, \lambda=1+2\sqrt{2} (transição caos-ordem em uma janela de período 3) e \lambda=4.0 (ergodicidade comprovada para esse valor de parâmetro). Considere os pontos do mapa logístico com precisão de 5 casas decimais, faça bom uso das escalas linear e logarítmica para que se possa observar a diferença entre cada resultado. (b) Calcule a função de correlação para os dígitos de pi e bases de nucleotídeos considerando apenas valores -1 e 1 para digitos abaixo e acima de 5, respectivamente, e para nucleotideos abaixo e acima de 2, respectivamente. === Lista 9: Automata celulares - 08/11/2011 === Para ler sobre automata celulares:\\ http://mathworld.wolfram.com/ElementaryCellularAutomaton.html\\ http://classes.yale.edu/fractals/CA/welcome.html\\ Desenvolva um código para receber o numero de uma dada regra de Wolfram como parâmetro de entrada e retornar a evolução temporal de uma condição inicial com 1 sítio central 1 os outros 0 por T passos de tempo em um sistema com N sítios.\\ (1) Mostre graficamente as evoluções temporais para as regras 30, 60, 90, 102, 126 e 150 com N=T=200. Utilize como condição inicial um único sítio (central) com estado s=1 e todos os outros com s=0 e implemente (a) condições de contorno periódicas (s[N]=s[0]) e (b) livres (s[N]=s[0]=0). \\ (2) Dada uma regra de automata (repita para as regras 30, 90, 102 e 126), condições de contorno nula e uma condição inicial com 1 sítio central com s=1, escolha um sítio i e acompanhe sua evolução temporal s(i,t) com t=0...T \\ a) Calcule o número médio de bits 1 ao longo da evolução com N=100000 e T=100, 1000, 5000, 10000, 50000 e 100000.\\ b) Calcule a função de correlação C(\tau) desta série de bits\\ (3) Agrupe os bits na evolução temporal do sitio i em blocos de 8 bits. Transforme o bloco de bits em um inteiro de 8 bits **com sinal**. \\ a) Faça um gráfico com a evolução temporal dos números de 8 bits. \\ b) Calcule o valor médio do número obtido. \\ c) Calcule a função de correlação para a série gerada com tal procedimento.\\ -- Para lembrar da representação de inteiros com sinal leia:\\ [[http://www.cs.cornell.edu/~tomf/notes/cps104/twoscomp.html|]]\\ [[http://en.wikipedia.org/wiki/Two's_complement|]] === Lista 10: Números pseudo-aleatórios e Monte-Carlo- 17/11/2011 === Números gerados com o computador, uma [[http://plato.stanford.edu/entries/turing-machine/|máquina de Turing determinística]], não pode gerar aleatoriedade com seus elementos lógicos. Podemos, no máximo, gerar sequências de números que passem por determinados [[http://www.stat.fsu.edu/pub/diehard/|testes estatísticos de aleatoriedade]]. A estes números chamamos //pseudo-aleatórios// e sua obtenção depende crucialmente do aumento de entropia gerado pela perda de bits (informação) nas operações lógicas - para os mais interessados: http://www.springerlink.com/content/jn7x3365386phn46/ Existem geradores de números verdadeiramente aleatórios, baseados na estocasticidade de processos físicos como emissão de fótons em processos radiativos ou espalhamento em divisores de feixes (veja [[http://qrbg.irb.hr/|aqui]] um exemplo). Esses geradores são usados para operações críticas, como criptografia e bingo. Uma boa fonte de informações a respeito de números pseudo-aleatórios é a biblioteca GSL, que possui, dentre diversas outras aplicações científicas como cálculo de autovalores e solução de sistemas de equações lineares, rotinas que geram números pseudo-aleatórios com diversas receitas. http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html - Usando o gerador congruencial x=(x*A)mod(B) gere 1000 pontos e visualize seu resultado em um espaço 2D considerando pares de numeros aleatórios consecutivos representando pontos no espaço. Utilize A=16807 e 65539 e B=2147483648. Estes números podem ser gerados sem a operação mod com números inteiros de 32 bits **com sinal.** Visualize seus resultados normalizando o número aleatório para gerar pontos no alcance [-1...1] - Usando o mesmo gerador com A=16807 e A=65539, crie pontos em um espaço 3D com cada ponto no espaço obtido a partir de uma trinca de números aleatórios gerados consecutivamente. Visualize seus resultados e, girando a figura no gnuplot, tente visualizar a regularidade nos pontos gerados com 65539. Esses pontos, de acordo com a [[http://www.pnas.org/content/61/1/25.citation|teoria]], caem em 15 planos! (Ver mais [[http://portal.acm.org/citation.cfm?id=363827|aqui]]. - Em 2D, visualize x versus 9x-6y+z onde x,y,z são a trinca de números aleatórios correspondente a um ponto no espaço 3D. Veja como o gerador 65539 é ruim. - Gere 10^4 pares de pontos (x,y) aleatoriamente distribuídos em 2D (com x e y entre -1 e 1). Calcule a fração de números gerados que satisfazem x^2+y^2 <=1 . Essa razão deve ser aproximadamente a razão entre a área de um círculo de raio unitário e um quadrado de lado 2. Esta razão vale \pi/4. Multiplique sua fração de números por 4 e voce terá calculado pi pelo Método de Monte Carlo. Faça um gráfico mostrando a convergência do valor estimado com o método de Monte Carlo para o valor aproximado de 3.14159265358979323846 (por exemplo, mostre no gráfico a diferença entre o valor aproximado e o valor estimado) utilizando o gerador com A=16807 e compare o resultado obtido com o gerador A=65539. Quantos números aleatórios devemos gerar para se obter 10\% de precisão na sua estimativa com A=16807? - Ruína do jogador: Considere dois jogadores, o primeiro com n1 moedas e o segundo com n2 moedas em um jogo de cara ou coroa. Cada vez que um jogador perde, ele entrega uma moeda para o outro jogador e o jogo se repete até um dos jogadores perder todo o dinheiro. Realize esse jogo um milhão de vezes e calcule a fração de vezes que cada jogador perdeu todo o dinheiro. Esse valor deve ser interpretado como a probabilidade P_i de cada jogador ir à falência. Compare seu resultado com os valores teóricos para as P_1=\frac{n_2}{n_1+n_2} e P_2=\frac{n_1}{n_1+n_2}. Para ler mais sobre o assunto [[http://mathworld.wolfram.com/GamblersRuin.html|clique aqui]]